#### installing packages if not installed ####

list.of.packages <- c('tidyverse',
                      'ggplot2',
                      'lubridate',
                      'estimatr',
                      'ggthemes',
                      'readr',
                      'readxl',
                      'rdrobust',
                      'gridExtra')
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

#### attaching libraries ####

suppressPackageStartupMessages(
  
  {
    
    library(tidyverse)
    library(meta)
    library(ggplot2)
    library(lubridate)
    library(estimatr)    
    library(ggthemes)
    library(readr)
    library(readxl)
    library(rdrobust)
    library(gridExtra)
    
  }
  
)

options(scipen=0)

#### loading in policing data for descriptive plots #### 

# Seattle
# 1) Terry Stops
# 2) Officer-initiated calls 

# loading terry stop data 

load("data/terry.RData")

terry = terry %>% 
  mutate(blm = ifelse(date >= as.Date("2020-05-29"), 1, 0)) %>% 
  mutate(initiate_off = ifelse(call_type == "-" | call_type == "ONVIEW" | 
                                 call_type == "SCHEDULED EVENT (RECURRING)" | 
                                 call_type == "ALARM CALL (NOT POLICE ALARM)", 1, 0),
         initiate_civ = 
           ifelse(call_type == "911" | 
                    call_type == "TELEPHONE OTHER, NOT 911" | 
                    call_type == "TEXT MESSAGE", 1, 0))

# loading 911 call data

load(file = "data/calldata.RData")

call = 
  call %>% 
  mutate(blm = ifelse(date >= as.Date("2020-05-29"), 1, 0)) %>% 
  mutate(initiate_off = ifelse(call_type == "-" | call_type == "ONVIEW" | 
                                 call_type == "SCHEDULED EVENT (RECURRING)" | 
                                 call_type == "ALARM CALL (NOT POLICE ALARM)" | 
                                 call_type == "PROACTIVE (OFFICER INITIATED)" | 
                                 call_type == "POLICE (VARDA) ALARM", 1, 0),
         initiate_civ = 
           ifelse(call_type == "911" | 
                    call_type == "TELEPHONE OTHER, NOT 911" | 
                    call_type == "TEXT MESSAGE" | 
                    call_type == "HISTORY CALL (RETRO)" | 
                    call_type == "IN PERSON COMPLAINT", 1, 0))

# gen time data

call$time_arrived2 = parse_date_time(call$time_arrived, "%m/%d/%Y %H:%M:%S %p")
call$time_queued2 = parse_date_time(call$time_queued, "%m/%d/%Y %H:%M:%S %p")

call$time_arrived2 = 
  as.POSIXct(call$time_arrived, format="%m/%d/%Y %H:%M:%S %p", tz="UTC")
call$time_queued2 = 
  as.POSIXct(call$time_queued, format="%m/%d/%Y %H:%M:%S %p", tz="UTC")
call$response_time = (call$time_arrived2 - call$time_queued2) / 60
the_index = call$response_time < 0
call = call %>% filter(!the_index)

# Philly 
# 1) Pedestrian Stops 
# 2) Traffic Stops 

# loading in stop data 

load(file = "data/ph_stops_ped.RData")
load(file = "data/ph_stops_veh.RData")

# Austin
# 1) Traffic Stops 

# loading in traffic stop data 

load(file = "data/rpf.RData")

# Los Angeles
# 1) Pedestrian stops
# 2) Traffic Stops
# 3) Officer-Initiated Calls

# stops 

load("data/stopdata_la.RData")

# calls

load("data/calls_la.RData")


#### rdit coefficients #### 

# Austin (Traffic Stops)

rpf_ts = rpf %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

rpf_ts$blmrv = 
  rpf_ts$date - min(rpf_ts$date[rpf_ts$blm == 1])

rpf_ts$count_std = 
  (rpf_ts$count - mean(rpf_ts$count, na.rm = TRUE)) / 
  sd(rpf_ts$count, na.rm = TRUE)

out1 = with(rpf_ts, rdrobust(x = blmrv, y = count_std,
                      kernel = "uni", p = 1, all = TRUE))

out_meantest = with(rpf_ts, rdrobust(x = blmrv, y = count,
                                           kernel = "uni", p = 1, all = TRUE))

abs(out_meantest$coef[1] / mean(rpf_ts$count[rpf_ts$blm == 0]))

# Los Angeles (Pedestrian Stops) 

stops_ped_ts = stopdat1 %>% 
  filter(traffic_stop == 0) %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

stops_ped_ts$blmrv = 
  stops_ped_ts$date - min(stops_ped_ts$date[stops_ped_ts$blm == 1])

stops_ped_ts$count_std = 
  (stops_ped_ts$count - mean(stops_ped_ts$count, na.rm = TRUE)) / 
  sd(stops_ped_ts$count, na.rm = TRUE)

out2 = with(stops_ped_ts, rdrobust(x = blmrv, y = count_std,
                             kernel = "uni", p = 1, all = TRUE))

out_meantest = with(stops_ped_ts, rdrobust(x = blmrv, y = count,
                                   kernel = "uni", p = 1, all = TRUE))

abs(out_meantest$coef[1] / mean(stops_ped_ts$count[stops_ped_ts$blm == 0]))
summary(out_meantest)

# Los Angeles (Traffic Stops) 

stops_veh_ts = stopdat1 %>% 
  filter(traffic_stop == 1) %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) 

stops_veh_ts$blmrv = 
  stops_veh_ts$date - min(stops_veh_ts$date[stops_veh_ts$blm == 1])

stops_veh_ts$count_std = 
  (stops_veh_ts$count - mean(stops_veh_ts$count, na.rm = TRUE)) / 
  sd(stops_veh_ts$count, na.rm = TRUE)

out3 = 
  with(stops_veh_ts, rdrobust(x = blmrv, y = count_std,
                              kernel = "uni", p = 1, all = TRUE))

out_meantest = with(stops_veh_ts, rdrobust(x = blmrv, y = count,
                                           kernel = "uni", p = 1, all = TRUE))

abs(out_meantest$coef[1] / mean(stops_veh_ts$count[stops_veh_ts$blm == 0]))

# Los Angeles (Officer-Initiated Calls) 

calls_off_ts = lacallsf %>% 
  mutate(count = 1) %>%
  filter(off_call == 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) 

calls_off_ts$blmrv = 
  calls_off_ts$date - min(calls_off_ts$date[calls_off_ts$blm == 1])

calls_off_ts$count_std = 
  (calls_off_ts$count - mean(calls_off_ts$count, na.rm = TRUE)) / 
  sd(calls_off_ts$count, na.rm = TRUE)

out4 = 
  with(calls_off_ts, rdrobust(x = blmrv, y = count_std,
                              kernel = "uni", p = 1, all = TRUE))

# Philly (Pedestrian Stops)

stops_ts = stops_ped %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) 

stops_ts$blmrv = 
  stops_ts$date - min(stops_ts$date[stops_ts$blm == 1])

stops_ts = 
  stops_ts %>% filter(date >= as.Date("2018-01-01"))

stops_ts$count_std = 
  (stops_ts$count - mean(stops_ts$count, na.rm = TRUE)) / 
  sd(stops_ts$count, na.rm = TRUE)

out5 = 
  with(stops_ts, rdrobust(x = blmrv, y = count_std,
                              kernel = "uni", p = 1, all = TRUE))

out_meantest = with(stops_ts, rdrobust(x = blmrv, y = count,
                                           kernel = "uni", p = 1, all = TRUE))

abs(out_meantest$coef[1] / mean(stops_ts$count[stops_ts$blm == 0]))

# Philly (Traffic Stops)

stops_ts = stops_veh %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) 

stops_ts$blmrv = 
  stops_ts$date - min(stops_ts$date[stops_ts$blm == 1])

stops_ts = 
  stops_ts %>% filter(date >= as.Date("2018-01-01"))

stops_ts$count_std = 
  (stops_ts$count - mean(stops_ts$count, na.rm = TRUE)) / 
  sd(stops_ts$count, na.rm = TRUE)

out6 = 
  with(stops_ts, rdrobust(x = blmrv, y = count_std,
                          kernel = "uni", p = 1, all = TRUE))

out_meantest = with(stops_ts, rdrobust(x = blmrv, y = count,
                                       kernel = "uni", p = 1, all = TRUE))

abs(out_meantest$coef[1] / mean(stops_ts$count[stops_ts$blm == 0]))

# seattle (terry stops)

terry_desc = terry %>% 
  filter(initiate_off == 1) %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm)) 

terry_desc$blmrv = 
  terry_desc$date - min(terry_desc$date[terry_desc$blm == 1])

terry_desc$count_std = 
  (terry_desc$count - mean(terry_desc$count, na.rm = TRUE)) / 
  sd(terry_desc$count, na.rm = TRUE)

out7 = 
  with(terry_desc, rdrobust(x = blmrv, y = count_std,
                          kernel = "uni", p = 1, all = TRUE))


out_meantest = with(terry_desc, rdrobust(x = blmrv, y = count,
                                       kernel = "uni", p = 1, all = TRUE))

abs(out_meantest$coef[1] / mean(terry_desc$count[terry_desc$blm == 0]))


# b) officer initiated calls 

call_desc = call %>% 
  filter(initiate_off == 1) %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count),
            blm = mean(blm))

call_desc$blmrv = 
  call_desc$date - min(call_desc$date[call_desc$blm == 1])

call_desc$count_std = 
  (call_desc$count - mean(call_desc$count, na.rm = TRUE)) / 
  sd(call_desc$count, na.rm = TRUE)

out8 = 
  with(call_desc, rdrobust(x = blmrv, y = count_std,
                            kernel = "uni", p = 1, all = TRUE))

# setting up figure

df_coef_plot = 
  data.frame(
    
    est = c(out1$coef[1], out2$coef[1], out3$coef[1], out4$coef[1],
            out5$coef[1], out6$coef[1], out7$coef[1], out8$coef[1]),
    
    se = c(out1$se[3], out2$se[3], out3$se[3], out4$se[3],
           out5$se[3], out6$se[3], out7$se[3], out8$se[3]),
    
    pv = c(out1$pv[3], out2$pv[3], out3$pv[3], out4$pv[3],
           out5$pv[3], out6$pv[3], out7$pv[3], out8$pv[3]),
    
    city = c("AUS", "LA", "LA", "LA", "PHL", "PHL", "SEA", "SEA"),
    
    outcome = factor(c("Traffic Stops", "Pedestrian Stops", "Traffic Stops", "Officer Calls",
                "Pedestrian Stops", "Traffic Stops", "Terry Stops", "Officer Calls"),
                levels = c("Traffic Stops", "Pedestrian Stops", "Officer Calls", "Terry Stops"))
    
  )


df_coef_plot_meta = 
  df_coef_plot %>% 
  filter(outcome != "Officer Calls") %>% 
  group_by(city) %>% 
  summarize(est = mean(est),
            se = mean(se))

meta_res = 
  metagen(TE = df_coef_plot_meta$est,
          seTE = df_coef_plot_meta$se, 
          sm = "SMD",
          method.random.ci = 'hk')


df_coef_plot_plot = df_coef_plot %>% 
  filter(outcome != "Officer Calls") %>% 
  bind_rows(., 
            data.frame(
              est = meta_res$TE.random,
              se = meta_res$seTE.random,
              pv = meta_res$pval.random,
              city = "Meta-Analysis",
              outcome = "Meta-Analysis"
            )
  ) %>% 
  mutate(city = factor(city,
                       levels = c("AUS","LA", "PHL", "SEA", "Meta-Analysis")),
         outcome = factor(outcome,
                          levels = c("Traffic Stops", "Pedestrian Stops",
                                     "Terry Stops", "Meta-Analysis")))

df_coef_plot_plot %>% 
  ggplot() + 
  geom_hline(yintercept = 0) + 
  geom_point(aes(x = city, y = est, shape = outcome),
             position = position_dodge(.35)) +
  geom_errorbar(aes(x = city, 
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    shape = outcome),
                width = 0,
                position = position_dodge(.35),
                size = .2) + 
  geom_errorbar(aes(x = city, 
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    shape = outcome),
                width = 0,
                position = position_dodge(.35),
                size = .4) + 
  labs(x = "City", y = "RDiT Coefficient\n(BLM Protest)",
       shape = "Outcome") + 
  theme_tufte(base_size = 9)

# analysis 

ggsave(plot = last_plot(), width = 6, height = 2, filename = "pics/figure2res.jpeg",
       dpi = 600)

#### rdit coefficients: placebos before and after #### 

# setting up the data 

# Austin (Traffic Stops)

rpf_ts = rpf %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

rpf_ts$blmrv = 
  rpf_ts$date - min(rpf_ts$date[rpf_ts$blm == 1])

rpf_ts$count_std = 
  (rpf_ts$count - mean(rpf_ts$count, na.rm = TRUE)) / 
  sd(rpf_ts$count, na.rm = TRUE)

# Los Angeles (Pedestrian Stops) 

stops_ped_ts = stopdat1 %>% 
  filter(traffic_stop == 0) %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

stops_ped_ts$blmrv = 
  stops_ped_ts$date - min(stops_ped_ts$date[stops_ped_ts$blm == 1])

stops_ped_ts$count_std = 
  (stops_ped_ts$count - mean(stops_ped_ts$count, na.rm = TRUE)) / 
  sd(stops_ped_ts$count, na.rm = TRUE)

# Los Angeles (Traffic Stops) 

stops_veh_ts = stopdat1 %>% 
  filter(traffic_stop == 1) %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) 

stops_veh_ts$blmrv = 
  stops_veh_ts$date - min(stops_veh_ts$date[stops_veh_ts$blm == 1])

stops_veh_ts$count_std = 
  (stops_veh_ts$count - mean(stops_veh_ts$count, na.rm = TRUE)) / 
  sd(stops_veh_ts$count, na.rm = TRUE)

out3 = 
  with(stops_veh_ts, rdrobust(x = blmrv, y = count_std,
                              kernel = "uni", p = 1, all = TRUE))

# Philly (Pedestrian Stops)

stops_ts_ped = stops_ped %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) 

stops_ts_ped$blmrv = 
  stops_ts_ped$date - min(stops_ts_ped$date[stops_ts_ped$blm == 1])

stops_ts_ped = 
  stops_ts_ped %>% filter(date >= as.Date("2018-01-01"))

stops_ts_ped$count_std = 
  (stops_ts_ped$count - mean(stops_ts_ped$count, na.rm = TRUE)) / 
  sd(stops_ts_ped$count, na.rm = TRUE)

# Philly (Traffic Stops)

stops_ts_veh = stops_veh %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) 

stops_ts_veh$blmrv = 
  stops_ts_veh$date - min(stops_ts_veh$date[stops_ts_veh$blm == 1])

stops_ts_veh = 
  stops_ts_veh %>% filter(date >= as.Date("2018-01-01"))

stops_ts_veh$count_std = 
  (stops_ts_veh$count - mean(stops_ts_veh$count, na.rm = TRUE)) / 
  sd(stops_ts_veh$count, na.rm = TRUE)

# seattle (terry stops)

terry_desc = terry %>% 
  filter(initiate_off == 1) %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm)) 

terry_desc$blmrv = 
  terry_desc$date - min(terry_desc$date[terry_desc$blm == 1])

terry_desc$count_std = 
  (terry_desc$count - mean(terry_desc$count, na.rm = TRUE)) / 
  sd(terry_desc$count, na.rm = TRUE)

terry_desc$blmrv

# generating loop for placebos

list_to_loop_df = 
  list(rpf_ts, stops_ped_ts, stops_veh_ts, stops_ts_ped, stops_ts_veh, terry_desc)
city_label = c("AUS", "LA", "LA", "PHL", "PHL", "SEA")
stop_label = c("Vehicle", "Pedestrian", "Vehicle", "Pedestrian", "Vehicle", "Terry")
list_to_loop_df_out = as.list(rep(NA, length(list_to_loop_df)))
seq_2_l = seq(from = -100, to = 100, by = 20)

for (i in 1:length(list_to_loop_df)) {
  
  print(paste0("I iteration ", i))
  
  df_seq_out = 
    matrix(NA, nrow = length(seq_2_l), ncol = 6) %>% 
    as.data.frame %>% 
    `colnames<-` (c("est", "se", "pv", 'plcval', 'city', 'outcome'))
  
  for (k in 1:length(seq_2_l)) {
    
    print(paste0("K iteration ", k))
    
    out = with(list_to_loop_df[[i]],
               rdrobust(x = I(blmrv + seq_2_l[k]), y = count, p = 1,
                        kernel = 'uni'))
    
    df_seq_out[k, 1] = out$coef[3]
    df_seq_out[k, 2] = out$se[3]
    df_seq_out[k, 3] = out$pv[3]
    df_seq_out[k, 4] = seq_2_l[k]
    df_seq_out[k, 5] = city_label[i]
    df_seq_out[k, 6] = stop_label[i]
    
  }
  
  list_to_loop_df_out[[i]] = df_seq_out
  
}

list_to_loop_df_out %>% 
  do.call(rbind.data.frame, .) %>% 
  mutate(city_stop = paste0(city, ', ', outcome)) %>% 
  mutate(plcvalbin = ifelse(plcval == 0, 1, 0)) %>% 
  mutate(plcvalbin = factor(plcvalbin, levels = c(0, 1))) %>% 
  ggplot() + 
  geom_hline(yintercept = 0) + 
  geom_point(aes(x = plcval, y = est,
                 col = plcvalbin),
             size = 2) + 
  geom_errorbar(aes(x = plcval, 
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    col = plcvalbin),
                size = .6,
                width = 0) + 
  geom_errorbar(aes(x = plcval, 
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    col = plcvalbin),
                size = .8,
                width = 0) + 
  guides(col = 'none') + 
  labs(x = "Placebo Discontinuity (Days to BLM Protest)",
       y = "RDiT Coefficient",
       title = 'Placebo Discontinuities Before and After the BLM Protest (Stop Outcome)') + 
  scale_color_grey(start = .6, end = 0) + 
  facet_wrap(~city_stop, scales = 'free') + 
  theme_tufte()

ggsave(plot = last_plot(), filename = 'pics/plcdisc_stops.jpeg', dpi = 600,
       width = 8, height = 4.5)

#### MAKE REGRESSION TABLE FOR APPENDIX ####

library(xtable)

atx_traffic <- as.data.frame(cbind(out1$coef[1],out1$se[3],out1$pv[3],
                                   out1$N[1],(out1$bws[1]*2),out1$bws[1]))

la_ped <- as.data.frame(cbind(out2$coef[1],out2$se[3],out2$pv[3],
                              out2$N[1],(out2$bws[1]*2),out2$bws[1]))

la_traffic <- as.data.frame(cbind(out3$coef[1],out3$se[3],out3$pv[3],
                              out3$N[1],(out3$bws[1]*2),out3$bws[1]))

ph_ped <- as.data.frame(cbind(out5$coef[1],out5$se[3],out5$pv[3],
                                  out5$N[1],(out5$bws[1]*2),out5$bws[1]))

ph_traffic <- as.data.frame(cbind(out6$coef[1],out6$se[3],out6$pv[3],
                                  out6$N[1],(out6$bws[1]*2),out6$bws[1]))

sea_terry <- as.data.frame(cbind(out7$coef[1],out7$se[3],out7$pv[3],
                                  out7$N[1],(out7$bws[1]*2),out7$bws[1]))

tab <- as.data.frame(rbind(atx_traffic,la_ped,la_traffic,
                           ph_ped,ph_traffic,sea_terry))

tab <- tab %>% mutate(V4 = as.character(V4))
tab <- tab %>% mutate_if(is.numeric,round, digits = 2)
tab <- tab %>% mutate(V3 = "0.000")

cities <- as.data.frame(rbind("Austin","LA","LA","Philly","Philly","Seattle"))
outcome <- as.data.frame(rbind("Traffic","Pedestrian","Traffic","Pedestrian","Traffic","Terry"))

tab <- as.data.frame(cbind(cities,outcome,tab))
colnames(tab) <- c("City","Stops","Coeff","SE","P-Val","N-Val","Effective N","Bandwidth Est.")

xt <- xtable(
  x = tab,
  align = "lllcccccc",
  type = "html",
  caption = "RDiT Coefficients Characterizing the Effect of BLM Protests on Policing Activities.",
  label = "tab:depol_reg_tab"
)


addtorow <- list()
addtorow$pos <- list(nrow(xt))
addtorow$command <- as.vector(c(paste("\\hline\n\\multicolumn{8}{l}{\\emph{All estimates are specified with a uniform kernel and polynomial degree equal to 1.}}\\\\\n",
                                      paste("\\multicolumn{8}{l}{\\emph{Standard errors are robust.}}\\\\\n",
                                      sep=""))))



print(xt,
      format.args = list(big.mark = ","),
      caption.placement = "top",
      include.rownames = FALSE,
      add.to.row = addtorow,
      file = "fig2_reg_tab.tex",
      hline.after = c(-1, 0, 6)
)

#### jacob blake effect #### 

# Austin (Traffic Stops)

rpf_ts = rpf %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

rpf_ts = rpf_ts %>% mutate(
  blake = ifelse(date >= as.Date("2020-08-23"), 1, 0)
) %>% 
  mutate(
    blakerv = date - as.Date('2020-08-23')
  )

rpf_ts$blmrv = 
  rpf_ts$date - min(rpf_ts$date[rpf_ts$blm == 1])

rpf_ts$count_std = 
  (rpf_ts$count - mean(rpf_ts$count, na.rm = TRUE)) / 
  sd(rpf_ts$count, na.rm = TRUE)

out1 = with(rpf_ts, rdrobust(x = blakerv, y = count_std,
                             kernel = "uni", p = 1, all = TRUE))

out_meantest = with(rpf_ts, rdrobust(x = blmrv, y = count,
                                     kernel = "uni", p = 1, all = TRUE))

abs(out_meantest$coef[1] / mean(rpf_ts$count[rpf_ts$blm == 0]))

# Los Angeles (Pedestrian Stops) 

stops_ped_ts = stopdat1 %>% 
  filter(traffic_stop == 0) %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE))

stops_ped_ts$blmrv = 
  stops_ped_ts$date - min(stops_ped_ts$date[stops_ped_ts$blm == 1])

stops_ped_ts$count_std = 
  (stops_ped_ts$count - mean(stops_ped_ts$count, na.rm = TRUE)) / 
  sd(stops_ped_ts$count, na.rm = TRUE)

stops_ped_ts = stops_ped_ts %>% mutate(
  blake = ifelse(date >= as.Date("2020-08-23"), 1, 0)
) %>% 
  mutate(
    blakerv = date - as.Date('2020-08-23')
  )


out2 = with(stops_ped_ts, rdrobust(x = blakerv, y = count_std,
                                   kernel = "uni", p = 1, all = TRUE))

out_meantest = with(stops_ped_ts, rdrobust(x = blakerv, y = count,
                                           kernel = "uni", p = 1, all = TRUE))

abs(out_meantest$coef[1] / mean(stops_ped_ts$count[stops_ped_ts$blm == 0]))
summary(out_meantest)

# Los Angeles (Traffic Stops) 

stops_veh_ts = stopdat1 %>% 
  filter(traffic_stop == 1) %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) 

stops_veh_ts$blmrv = 
  stops_veh_ts$date - min(stops_veh_ts$date[stops_veh_ts$blm == 1])

stops_veh_ts$count_std = 
  (stops_veh_ts$count - mean(stops_veh_ts$count, na.rm = TRUE)) / 
  sd(stops_veh_ts$count, na.rm = TRUE)

stops_veh_ts = stops_veh_ts %>% mutate(
  blake = ifelse(date >= as.Date("2020-08-23"), 1, 0)
) %>% 
  mutate(
    blakerv = date - as.Date('2020-08-23')
  )

out3 = 
  with(stops_veh_ts, rdrobust(x = blakerv, y = count_std,
                              kernel = "uni", p = 1, all = TRUE))

out_meantest = with(stops_veh_ts, rdrobust(x = blakerv, y = count,
                                           kernel = "uni", p = 1, all = TRUE))

abs(out_meantest$coef[1] / mean(stops_veh_ts$count[stops_veh_ts$blm == 0]))


# Philly (Pedestrian Stops)

stops_ts = stops_ped %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) 

stops_ts$blmrv = 
  stops_ts$date - min(stops_ts$date[stops_ts$blm == 1])

stops_ts = 
  stops_ts %>% filter(date >= as.Date("2018-01-01"))

stops_ts$count_std = 
  (stops_ts$count - mean(stops_ts$count, na.rm = TRUE)) / 
  sd(stops_ts$count, na.rm = TRUE)

stops_ts = stops_ts %>% mutate(
  blake = ifelse(date >= as.Date("2020-08-23"), 1, 0)
) %>% 
  mutate(
    blakerv = date - as.Date('2020-08-23')
  )

out5 = 
  with(stops_ts, rdrobust(x = blakerv, y = count_std,
                          kernel = "uni", p = 1, all = TRUE))

out_meantest = with(stops_ts, rdrobust(x = blmrv, y = count,
                                       kernel = "uni", p = 1, all = TRUE))

abs(out_meantest$coef[1] / mean(stops_ts$count[stops_ts$blm == 0]))

# Philly (Traffic Stops)

stops_ts = stops_veh %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm, na.rm = TRUE)) 

stops_ts$blmrv = 
  stops_ts$date - min(stops_ts$date[stops_ts$blm == 1])

stops_ts = 
  stops_ts %>% filter(date >= as.Date("2018-01-01"))

stops_ts$count_std = 
  (stops_ts$count - mean(stops_ts$count, na.rm = TRUE)) / 
  sd(stops_ts$count, na.rm = TRUE)

stops_ts = stops_ts %>% mutate(
  blake = ifelse(date >= as.Date("2020-08-23"), 1, 0)
) %>% 
  mutate(
    blakerv = date - as.Date('2020-08-23')
  )


out6 = 
  with(stops_ts, rdrobust(x = blakerv, y = count_std,
                          kernel = "uni", p = 1, all = TRUE))

out_meantest = with(stops_ts, rdrobust(x = blakerv, y = count,
                                       kernel = "uni", p = 1, all = TRUE))

abs(out_meantest$coef[1] / mean(stops_ts$count[stops_ts$blm == 0]))

# seattle (terry stops)

terry_desc = terry %>% 
  filter(initiate_off == 1) %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count, na.rm = TRUE),
            blm = mean(blm)) 

terry_desc$blmrv = 
  terry_desc$date - min(terry_desc$date[terry_desc$blm == 1])

terry_desc$count_std = 
  (terry_desc$count - mean(terry_desc$count, na.rm = TRUE)) / 
  sd(terry_desc$count, na.rm = TRUE)

terry_desc = terry_desc %>% mutate(
  blake = ifelse(date >= as.Date("2020-08-23"), 1, 0)
) %>% 
  mutate(
    blakerv = date - as.Date('2020-08-23')
  )

out7 = 
  with(terry_desc, rdrobust(x = blakerv, y = count_std,
                            kernel = "uni", p = 1, all = TRUE))

out_meantest = with(terry_desc, rdrobust(x = blakerv, y = count,
                                         kernel = "uni", p = 1, all = TRUE))


# b) officer initiated calls 

call_desc = call %>% 
  filter(initiate_off == 1) %>% 
  mutate(count = 1) %>% 
  group_by(date) %>% 
  summarize(count = sum(count),
            blm = mean(blm))

call_desc$blmrv = 
  call_desc$date - min(call_desc$date[call_desc$blm == 1])

call_desc$count_std = 
  (call_desc$count - mean(call_desc$count, na.rm = TRUE)) / 
  sd(call_desc$count, na.rm = TRUE)

out8 = 
  with(call_desc, rdrobust(x = blmrv, y = count_std,
                           kernel = "uni", p = 1, all = TRUE))

# setting up figure

df_coef_plot = 
  data.frame(
    
    est = c(out1$coef[1], out2$coef[1], out3$coef[1], out4$coef[1],
            out5$coef[1], out6$coef[1], out7$coef[1], out8$coef[1]),
    
    se = c(out1$se[3], out2$se[3], out3$se[3], out4$se[3],
           out5$se[3], out6$se[3], out7$se[3], out8$se[3]),
    
    pv = c(out1$pv[3], out2$pv[3], out3$pv[3], out4$pv[3],
           out5$pv[3], out6$pv[3], out7$pv[3], out8$pv[3]),
    
    city = c("AUS", "LA", "LA", "LA", "PHL", "PHL", "SEA", "SEA"),
    
    outcome = factor(c("Traffic Stops", "Pedestrian Stops", "Traffic Stops", "Officer Calls",
                       "Pedestrian Stops", "Traffic Stops", "Terry Stops", "Officer Calls"),
                     levels = c("Traffic Stops", "Pedestrian Stops", "Officer Calls", "Terry Stops"))
    
  )


df_coef_plot_meta = 
  df_coef_plot %>% 
  filter(outcome != "Officer Calls") %>% 
  group_by(city) %>% 
  summarize(est = mean(est),
            se = mean(se))

meta_res = 
  metagen(TE = df_coef_plot_meta$est,
          seTE = df_coef_plot_meta$se, 
          sm = "SMD",
          method.random.ci = 'hk')


df_coef_plot_plot = df_coef_plot %>% 
  filter(outcome != "Officer Calls") %>% 
  bind_rows(., 
            data.frame(
              est = meta_res$TE.random,
              se = meta_res$seTE.random,
              pv = meta_res$pval.random,
              city = "Meta-Analysis",
              outcome = "Meta-Analysis"
            )
  ) %>% 
  mutate(city = factor(city,
                       levels = c("AUS","LA", "PHL", "SEA", "Meta-Analysis")),
         outcome = factor(outcome,
                          levels = c("Traffic Stops", "Pedestrian Stops",
                                     "Terry Stops", "Meta-Analysis")))

df_coef_plot_plot %>% 
  ggplot() + 
  geom_hline(yintercept = 0) + 
  geom_point(aes(x = city, y = est, shape = outcome),
             position = position_dodge(.35)) +
  geom_errorbar(aes(x = city, 
                    ymin = est - 1.96 * se,
                    ymax = est + 1.96 * se,
                    shape = outcome),
                width = 0,
                position = position_dodge(.35),
                size = .2) + 
  geom_errorbar(aes(x = city, 
                    ymin = est - 1.645 * se,
                    ymax = est + 1.645 * se,
                    shape = outcome),
                width = 0,
                position = position_dodge(.35),
                size = .4) + 
  labs(x = "City", y = "RDiT Coefficient\n(BLM Protest)",
       shape = "Outcome") + 
  theme_tufte(base_size = 9)

ggsave(plot = last_plot(), width = 6, height = 2, filename = "pics/figure2resblake.jpeg",
       dpi = 600)
